   setwd("C:/Users/Jef/Desktop/Applied Statistics")
#   setwd("D:/Mijn Documenten/applied statistics - project")
library(qcc)

# opgave 3b   

   # we kiezen v=30 omdat dit een extreme waarde geeft!! --> uitleggen dat zit de maximale waarde geeft voor Hplus --> blablabla
v <- 30
   # door verschillende kritieke waarden te testen zoeken we die Hplus zo dat ARL is ongeveer 200
# Hplus <- 6.5*sqrt(1.4)
Hplus <- 7.60
RL <- array(dim=c(10000,1))
   # we willen 10000 observaties van de RL om de ARL te bepalen
for(i in 1:10000){
  j=0
  Tj=0
  while(Tj <= Hplus){
    j <- j + 1
    xj <- rnorm(1, mean=v, sd=sqrt(1.4))
    Wj <- -1/(2*(1.4^2))*(  (xj-max(xj,30))^2-(xj-min(xj,30))^2  )
    Tj <- Wj + max(0,Tj)
  }   
  # we willen voor alle i de j bijhouden:
  RL[i,1] <- j
}
ARL <- mean(RL)
hist(RL)
boxplot(RL)
ARL
  # opmerking: ARL is niet erg stabiel, dit komt doordat de RL waarden exponentieel verdeeld zijn (zie plaatje gemaakt door "hist(RL)")
  # dit beinvloedt dus ook erg de waarde voor je ARL !! 
  # (discussie over weggooien outliers ??? --> niet doen, wel opmerking over maken!)
  # en ook: wat voor effect heeft dit op vraag c?    (bij vraag 3c)
  
  # waarom staat in dictaat dat het numeriek moet --> het kan ook analitisch !! --> uitwerken!!


# opgave 3c
viscosity <- read.table("viscosity.txt")
Tvalues <- array(dim=c(170,1))
Tj <- 0

for(j in 1:170){   
  xj <- viscosity[j,1]  
  Wj <- -1/(2*(1.4^2))*(  (xj-max(xj,30))^2-(xj-min(xj,30))^2  )
  Tj <- Wj + max(0,Tj)
  Tvalues[j,1] <- Tj   
}
Tvalues
Tvalues < 7.65
# dus hij geeft een out-of-control melding vanaf de 94e observatie